
package org.zhongyin.coolingload;
import java.io.IOException;
import java.text.ParseException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import jxl.read.biff.BiffException;

import Jama.Matrix;


public class Regression {
	
	private static final double TUNING_CONSTANT = 4.685d;
	
	private Matrix coefficient;
	
	private double leastSquares;
	
	private double standardDeviation;
	
	private Matrix weight;
	
	//D = sqrt(eps(class(double))) from STATROBUSTFIT_Revised.m under matlab
	private static final double D = 1.490116119384766e-08;
	
	private static final int iterLimit = 50;
	
	public Regression(){
		coefficient = null;
		
		leastSquares = 0d;
		
		standardDeviation = 0d;
		
		weight = null;
		
	}
	
	
	
	/**
	 * For Linear Regression
	 * b = (XT*X)-1 * XT * Y
	 * @param array_X	X
	 * @param array_Y	Y
	 * @return	coefficient (b)
	 */
	public Matrix getCoefficient(double[][] array_X, double[] array_Y){
		
		if(array_X == null || array_Y == null){
			throw new NullPointerException();
		}
		Matrix X = new Matrix(array_X);
		Matrix Y = new Matrix(array_Y, array_Y.length);
		
		Matrix X_Transpose = X.transpose();
		// XT*X
		Matrix X_Transpose_Times_X = X_Transpose.times(X);
		//(XT*X)-1 Matrix inverse
		Matrix square_inverse = X_Transpose_Times_X.inverse();
		//(XT*X)-1 * XT
		Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
		
		Matrix result = inverseTimesTranspose.times(Y);
		
		
		coefficient = result;
		
		return result;
		
	}
	
	
	// for linear regression
	/**
	 * @param array_Y	Y
	 * @param array_X	X
	 * @param array_b	coefficient
	 * @return	leastsquares (e2)
	 */
	public double getLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b){
		if(array_Y == null || array_X == null || array_b == null){
			throw new NullPointerException();
		}
		
		Matrix X = new Matrix(array_X);
		Matrix b = new Matrix(array_b);
		Matrix Xb = X.times(b);
		
		
		double[][] array_Xb = Xb.getArray();
		
		double[][] e = new double[array_Y.length][1];
		
		
		double result = 0d;
		
		
		for(int i = 0; i < e.length; i++){
			
			e[i][0] = array_Y[i] - array_Xb[i][0];
			result += e[i][0] * e[i][0];
		}
		
		
		
		leastSquares = result;
		
		return result;
		
		
	}
	
	
	/**
	 * ri = (yi - ybari) / std
	 * @param array_Y 	Y
	 * @param array_X	X
	 * @param array_b	coefficient
	 * @param sDeviation 
	 * @return	Standardized residuals
	 */
	public Matrix getStdResidual(double[] array_Y, double[][] array_X, double[][] array_b, double sDeviation){
		if(array_Y == null || array_X == null || array_b == null){
			throw new NullPointerException();
		}
		
		Matrix X = new Matrix(array_X);
		Matrix b = new Matrix(array_b);
		Matrix Xb = X.times(b);
		
		double[][] y_bar = Xb.getArray();
		
		double[] r = new double[array_Y.length];
		
		for(int i = 0; i < r.length; i++){
			r[i] = (array_Y[i] - y_bar[i][0])/sDeviation;
		}
		
		Matrix stdResidual = new Matrix(r, array_Y.length);
		
		return stdResidual;
	}
	
	
	
	/**
	 * 
	 * Compute the adjusted residuals
	 * @param array_Y Y
	 * @param array_X X
	 * @param array_b coefficient
	 * @param array_h leverage
	 * @return
	 */
	public double[] getAdjustedResiduals(double[] array_Y, double[][] array_X, double[][] array_b, Matrix array_h){
		if(array_Y == null || array_X == null || array_b == null){
			throw new NullPointerException();
		}
		
		Matrix X = new Matrix(array_X);
		Matrix b = new Matrix(array_b);
		Matrix Xb = X.times(b);
		
		double[][] h = array_h.getArray();
		
		double[][] y_bar = Xb.getArray();
		
		double[] r = new double[array_Y.length];
		
		for(int i = 0; i < r.length; i++){
//			System.out.print(array_Y[i] - y_bar[i][0]+ ", ");
			r[i] = (array_Y[i] - y_bar[i][0])/Math.sqrt(1-h[i][i]);
//			r[i] = 1/Math.sqrt(1-h[i][i]);
		}
		
//		System.out.println();
		
		return r;
	
	}
	
	/**
	 * s is the robust variance given by MAD/0.6745 where MAD 
	 * is the median absolute deviation of the residuals
	 * 
	 * @param array_Radj  Radj
	 * @param array_X X
	 * @return
	 */
	public double getS(double[] array_Radj, double[][] array_X){
		if(array_Radj == null){
			return 0;
		}
		
		
		 double[] array_Radj_Clone = array_Radj.clone();
		
		 double median = 0d;
		 
	     int length = array_X[0].length; 
	     
//	     System.out.println ("length : " + length);
	     for(int i = 0; i < array_Radj.length; i++){
	    	 array_Radj_Clone[i] = Math.abs(array_Radj[i]);
	     }
		
		 Arrays.sort(array_Radj_Clone);
		 
		 
		 
		 array_Radj_Clone = Arrays.copyOfRange(array_Radj_Clone, length-1, array_Radj.length);
		 
//		 System.out.println("Sorted Radj: " + Arrays.toString(array_Radj_Clone));
		 
		 if(array_Radj.length % 2 == 0){
			 
			 
			 median = (array_Radj_Clone[array_Radj_Clone.length/2] + array_Radj_Clone[array_Radj_Clone.length/2 - 1])/2;
		 }
		 else{
			 median = array_Radj_Clone[array_Radj_Clone.length/2];
		 }
		 
//		 median = array_Radj[array_Radj.length/2];
		 
//		 System.out.println("Media :" + median);
		 double MAD = median;
		 
		 double S = MAD / 0.6745;
		 
		 return S;
		
		
	}
	
	/**
	 * 
	 * @param stdResidual ri
	 * @return	W[][] with Wii(ri)
	 */
	public Matrix getOriginalWeight(Matrix stdResidual){
		if(stdResidual == null){
			throw new NullPointerException();
		}
		
		double[][] r = stdResidual.getArray();
		
		double[][] W = new double[r.length][r.length];
		
		for(int i = 0; i < r.length; i++){
			
			if(Math.abs(r[i][0]) <= TUNING_CONSTANT){
				W[i][i] = (1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT))*(1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT));
			}
			else{
				W[i][i] = 0d;
			}
		}
		
		
		
		Matrix Weight = new Matrix(W);
		
		weight = Weight;
		
		return Weight;
		
		
	}
	
	public Matrix getRobustWeight(Matrix U){
		if(U == null){
			throw new NullPointerException();
		}
		
		double[][] u = U.getArray();
		
		double[][] W = new double[u.length][u.length];
		
		for(int i = 0; i < W.length; i++){
			
			if(Math.abs(u[i][0]) <= 1){
				W[i][i] = (1 - (u[i][0] * u[i][0]))*(1 - (u[i][0] * u[i][0]));
			}
			else{
				W[i][i] = 0d;
			}
		}
		
		
		
		Matrix Weight = new Matrix(W);
		
		weight = Weight;
		
		return Weight;
		
		
	}
	
	
	/**
	 * 
	 * The quantities hi are the leverages
	 * , the diagonal elements of X(Xt*X)-1 * Xt
	 * @param array_X X
	 * @return array_h
	 */
	public Matrix getLeverage(double[][] array_X){
		if(array_X == null){
			throw new NullPointerException();
		}
		
		Matrix X = new Matrix(array_X);
		
		// h = X(XTX)-1XT
		Matrix result = X.times((X.transpose().times(X)).inverse()).times(X.transpose());
		
		return result;
	}
	
	
	/**
	 * @param leastSquares e2
	 * @return e
	 */
	public double getStandardDeviation(double leastSquares){
		
		double result = Math.sqrt(leastSquares);
		
		standardDeviation = result;
		
		return result;
		
	}
	
	
	
	/**
	 * For Weighted Linear Regression
	 * b = (XT*WX)-1 * XT * W *Y
	 * @param array_X	X
	 * @param array_Y	Y
	 * @return	coefficient (b)
	 */
	public Matrix getWeightedCoefficient(double[] array_Y, double[][] array_X, Matrix weight){
		if(array_X == null || array_Y == null || weight == null){
			throw new NullPointerException();
		}
		Matrix X = new Matrix(array_X);
		Matrix Y = new Matrix(array_Y, array_Y.length);
		
		Matrix X_Transpose = X.transpose();
		
		//XT*W
		Matrix X_Transpose_Times_W = X_Transpose.times(weight);
		// XT*W*X
		Matrix X_Transpose_Times_W_Times_X = X_Transpose_Times_W.times(X);
		//(XT*WX)-1 Matrix inverse
		Matrix square_inverse = X_Transpose_Times_W_Times_X.inverse();
		//(XT*WX)-1 * XT
		Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
		//(XT*WX)-1 * XTW
		Matrix timesWeight = inverseTimesTranspose.times(weight);
		//(XT*WX)-1 * XTWY
		Matrix result = timesWeight.times(Y);
		
		
		coefficient = result;
		
		return result;
	}
	
	
	
	
	/**
	 * For Weighted Linear Regression
	 * @param array_Y	Y
	 * @param array_X	X
	 * @param array_b	coefficient
	 * @return	leastsquares (e2)
	 */
	public double getWeightedLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b, Matrix weight){
		if(array_Y == null || array_X == null || array_b == null || weight == null){
			throw new NullPointerException();
		}
		
		Matrix X = new Matrix(array_X);
		Matrix b = new Matrix(array_b);
		Matrix Xb = X.times(b);
		
		
		double[][] array_Xb = Xb.getArray();
		
		double[][] e = new double[array_Y.length][1];
		
		
		double result = 0d;
		
		
		for(int i = 0; i < e.length; i++){
			
			e[i][0] = array_Y[i] - array_Xb[i][0];
			//Wii*(Yi-Xib);
			result += weight.getArray()[i][i] * (e[i][0] * e[i][0]);
		}
		
		
		
		leastSquares = result;
		
		return result;
		
		
	}
	
	
	
	public Matrix getU(double[] array_Radj, double S){
		if(array_Radj == null){
			return null;
		}
		
		double[] u = new double[array_Radj.length];
		
		
		
		
		for(int i = 0; i < array_Radj.length; i++){
			u[i] = array_Radj[i]/(TUNING_CONSTANT*S);
		}
		
		Matrix U = new Matrix(u, u.length);
		
		return U;
		
		
	}
	/**
	 * Iteratively get the weight by Irls
	 * @param array_Y Y
	 * @param array_X X
	 * @param weight array of weight
	 * @return get more accurate weight
	 */
	public Matrix getIterWeight(double[] array_Y, double[][] array_X, Matrix weight){
		
//		Regression regression = new Regression();
		
		int length = weight.getArray().length;
		
		int iter = 0;
		
//		//Generate an identity matrix
//		weight = Matrix.identity(length, length);
		
		double[] weight_array = new double[length];
		
		for(int i = 0 ; i < length; i++){
			weight_array[i] = weight.get(i, i);
		}
		
//		System.out.println(iter + " weight: " + Arrays.toString(weight_array));
		// First coefficient (b)
		Matrix array_b = getWeightedCoefficient(array_Y, array_X, weight);
//		System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
		// First leverage
		Matrix leverage = getLeverage(array_X);
//		System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
		// First adjusted residuals
		double[] array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);
//		System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
		// First S
		double S = getS(array_Radj, array_X); 
//		System.out.println(iter + " S: " + S);
		// First U
		Matrix U = getU(array_Radj, S);
//		System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
		
		
		
		
		
//		System.out.println();
		
		
		
		while(any(array_b.getArray()) && iter < iterLimit){
			
			iter++;
			
			weight = getRobustWeight(U);
			
			for(int i = 0 ; i < length; i++){
				weight_array[i] = weight.get(i, i);
			}
			
//			System.out.println(iter + " weight: " + Arrays.toString(weight_array));
			
			array_b = getWeightedCoefficient(array_Y, array_X, weight);
//			System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
			// First leverage
			leverage = getLeverage(array_X);
//			System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
			// First adjusted residuals
			array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);
//			System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
			// First S
			S = getS(array_Radj, array_X); 
//			System.out.println(iter + " S: " + S);
			// First U
			U = getU(array_Radj, S);
//			System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
			
			
			
			
			
//			System.out.println();
			
			
		}
		
		
		double[][] diagonal_weight = new double[weight_array.length][weight_array.length];
		
		for(int i = 0 ; i < weight_array.length; i++){
			diagonal_weight[i][i] = weight_array[i];
		}
		Matrix result = new Matrix(diagonal_weight);
		
		return result;
	}
	
	/**
	 * 
	 * True if any element of a vector is
	 * logical 1 (TRUE).
	 * @param array
	 * @return
	 */
	public boolean any(double[][] array){
		if(array == null){
			return false;
		}
		
		for(int i = 0; i < array.length; i++){
			if(Math.abs(array[i][0])  <= D * Math.max(Math.abs(array[i][0]), 0)){
				return false;
			}
		}
		
		return true;
	}
	
	
	/**
	 * get forgetting factor based on lumda from 0.970 to  1.000
	 * @param weight
	 * @return a list of forgetting factor
	 */
	public List<double[][]> getWeightPrimePrime(Matrix weight){
		if(weight == null){
			throw new NullPointerException();
		}
		
		double k = 0.971d;
		
		List<double[][]> list = new ArrayList<double[][]>();
		
		int length = weight.getArray().length;
		
		
		while( k <= 0.971){
			Matrix lumda = Matrix.identity(length, length);
			double[][] array_lumda = lumda.getArray();
			
			System.out.print("w: ");
			for(int i = 0; i < length; i++){
				array_lumda[i][i] = lumda.get(i, i) * Math.pow(k, length-1 - i) * weight.get(i, 0);
				System.out.print(array_lumda[i][i]+ ", ");
			}
			

			list.add(array_lumda);
			k += 0.001;
		}
		
		return list;
	}
	
	public double getSecondStageWeightedLeastSquare(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
		if(array_Y == null || array_X == null || weightPrimePrime == null){
			throw new NullPointerException();
		}
		
		double minE = Double.MAX_VALUE;
		double minLeastSquare = Double.MAX_VALUE;
		
		int counter = 0;
		for(double[][] array_weight: weightPrimePrime){
			
			if(counter == 29){
				System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
			}
			
//			System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
			
			Matrix weight = new Matrix(array_weight);
			Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
			
			double[][] array_b = coefficient.getArray();
			
			System.out.println(counter+ "b:" + Arrays.deepToString(array_b));
			
			double weightedLeastSquare = getWeightedLeastSquares(array_Y, array_X, array_b, weight);
			
			double e = Math.sqrt(weightedLeastSquare/array_Y.length);
			
			System.out.println(counter+ "e2: " + weightedLeastSquare);
			
			System.out.println(counter+ "e: " + e);
			
			minE = Math.min(minE, e);
			
			minLeastSquare = Math.min(minLeastSquare, weightedLeastSquare);
			
			counter++;
			
		}
		
		return minE;
	}
	
	
	public List<Matrix> getFFCofficient(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
		if(array_Y == null || array_X == null || weightPrimePrime == null){
			throw new NullPointerException();
		}
		
		double lumbda = 0.971d;
		
		
		
		List<Matrix> result = new ArrayList<Matrix>(); 
		for(double[][] array_weight: weightPrimePrime){
			
			
			

			Matrix weight = new Matrix(array_weight);
			Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
			
			System.out.println();
			
			System.out.println(lumbda + " b: " + Arrays.deepToString(coefficient.getArray()));
			
			result.add(coefficient);
			
			
			lumbda += 0.001;
			
		}
		
		return result;
		
		
	}
	
	public static void main(String args[]) throws IOException, ParseException, BiffException{
////		double[][] array_X = {{1,2},{1,4},{1,7}};
////		double[] array_Y = {1,3,5};
////		double[] array_Weight = {1,1,1};
//////		
////		int length = 3;
//		
//		ReadExcel Original = new ReadExcel();
//		Original.setInputFile("c:/CoolingLoadProject/Original data.xls");
//		
//	    double[][] original_data = Original.read();
//	    
//	    int length = original_data.length;
//	    double[] array_Y = new double[original_data.length];
//	    double[][] array_X = new double[original_data.length][17];
//	    
//	    
//	    
//	    for(int i = 0; i < original_data.length; i++){
//	    	array_Y[i] = original_data[i][0];
//	    	for(int j = 0; j < 11; j++){	
//	    		array_X[i][j] = original_data[i][j+1];
//	    	}
//	    }
//	    
//	    for(int i = 0; i < original_data.length; i++){
//	    		//t*t
//	    		array_X[i][11] = original_data[i][6] * original_data[i][6];
//	    		//t24*t24
//	    		array_X[i][12] = original_data[i][7] * original_data[i][7];
//	    		//t168*t168
//	    		array_X[i][13] = original_data[i][8] * original_data[i][8];
//	    		//t * RHt
//	    		array_X[i][14] = original_data[i][6] * original_data[i][9];
//	    		//t24 *RHt24
//	    		array_X[i][15] = original_data[i][7] * original_data[i][10];
//	    		//t168 * RHt168
//	    		array_X[i][16] = original_data[i][8] * original_data[i][11];
//	    	
//	    }
//	    
//	    
//	    
////	    Matrix Y = new Matrix(array_Y, array_Y.length);
////	    
////	    double[] test_Y = Y.getRowPackedCopy();
////	    
////	    System.out.println(Arrays.toString(test_Y));
//	    
//	    
//	    
//	    Utility.printArray(array_X, 0, length);
	    
//	    System.out.println(Arrays.toString(array_Y));
		
		ARXWeekDays test = new ARXWeekDays();
		
		test.setARXmodel();
//		
		Utility.printArraySingleLine(test.getARXmodel_X());
		
		System.out.println(Arrays.toString(test.getARXmodel_Y()));
		
		double[][] array_X = test.getARXmodel_X();
		
		double[] array_Y = test.getARXmodel_Y();
		
		int length = array_Y.length;
	    
	    Matrix weight = Matrix.identity(length, length);
		
		Regression regression = new Regression();
		
		weight = regression.getIterWeight(array_Y, array_X, weight);
		
		System.out.println("IRLS Weight: " + Arrays.deepToString(weight.getArray()));
		
		Matrix cofficient = regression.getWeightedCoefficient(array_Y, array_X, weight);
		
		System.out.println("IRLS b: " + Arrays.deepToString(cofficient.getArray()));
		
		
//		List<double[][]> weightPrimePrime = new ArrayList<double[][]>();
		
//		weightPrimePrime = regression.getWeightPrimePrime(weight);
		
//		double minLeastSquare = regression.getSecondStageWeightedLeastSquare(array_Y, array_X, weightPrimePrime);
		
//		System.out.println("Min Least Square: " + minLeastSquare);
		
//		regression.getFFCofficient(array_Y, array_X, weightPrimePrime);
		
	}
	
	
}
